International Journal for Numerical Methods in Biomedical Engineering

Thin structures such as heart valves and aortic dissection flaps interact dynamically with blood flow in human vessels. Their flexibility and capacity for large deformations generate complex, highly transient hemodynamic patterns over the cardiac cycle. Accurately resolving these interactions remains challenging for conventional boundary-fitted fluid-structure interaction approaches. We present an immersed boundary method for simulating thin structures in incompressible flow on unstructured grids. The method couples a stabilized finite element fluid solver with a nonlinear, rotation-free shell formulation through a direct forcing immersed boundary approach. The framework supports both weak (explicit) and strong (implicit) time-coupling strategies, enabling stable simulations over a wide range of solid-to-fluid density ratios. Hydrodynamic forces acting on thin structures are computed from fluid solutions sampled on both sides of the structure, allowing accurate force reconstruction for zero-thickness shells. To our knowledge, this is the first immersed boundary formulation that couples an unstructured finite element fluid solver with a two-dimensional, rotation-free shell model to simulate interactions between thin structures and incompressible flow. Fluid-structure coupling is achieved using predefined finite element shape functions, which provide consistent projection between Eulerian and Lagrangian fields without additional interpolation procedures. The framework is validated using three-dimensional benchmark problems involving thin structures. Then, valve-like model is used to compare strong and weak coupling strategies. Finally, the method is applied to an idealized type-B aortic dissection model. The proposed approach is implemented within the open-source software CRIMSON, a finite element platform for cardiovascular simulation.

The mechanical behavior of right ventricular (RV) myocardium is governed by its anisotropic microstructure, yet constitutive models that account for fiber dispersion and enable reliable parameter identification remain limited. In this study, we propose a physics-embedded constitutive neural network framework for automated discovery of strain energy functions and microstructural parameters from experimental data. The model is formulated within an incompressible, orthotropic hyperelastic setting using invariant-based representations. Fiber, sheet, and normal directions are incorporated through a rotated structural basis, and dispersion effects are modeled using a generalized structure tensor approach. The framework is trained on multi-axial mechanical data from ovine RV myocardium, including uniaxial tension-compression and simple shear tests. We investigate two training scenarios: (i) full datasets containing both tensile and compressive regimes and (ii) datasets restricted to tensile loading. In both cases, the model accurately reproduces the measured stress-strain responses and identifies sparse, interpretable constitutive models which involve isotropic, anisotropic, and coupling invariants. However, the identifiability of microstructural parameters strongly depends on the available loading conditions. While tensile-only data yield higher predictive accuracy, they result in non-unique or biased estimates of fiber dispersion. In contrast, inclusion of compressive data enables consistent identification of dispersion parameters by separating fiber and matrix contributions. These results highlight the importance of multi-axial loading data for robust parameter identification and demonstrate the capability of constitutive neural network-based approaches for data-driven modeling of anisotropic soft tissues.

Aortic stiffness is both heterogenous and anisotropic. Current non-invasive methods to estimate aortic stiffness are limited to characterizing the aortic tissue as isotropic due to the lack the techniques required to extract multi-axial strain from 3D dynamic images. Vascular deformation mapping (VDM) is a nonrigid image registration technique which has thus far been applied to map aortic growth using longitudinal imaging. In this study, we propose to use VDM to assess 3D aortic deformation by mapping diastolic and systolic images. During image registration process, penalty parameters are employed to fine-tune image alignment and penalize non-physiological deformations. These penalty parameters must be calibrated to ensure that VDM successfully reproduces multi-axial aortic motion patterns in health and disease. In this paper, we developed a calibration pipeline for these parameters using synthetic data. A rotation-free shell model was used to generate physics-based synthetic data on aortic motion incorporating patient-specific geometries, root motion, and blood pressure from a cohort of 14 subjects (healthy, Marfans syndrome and thoracic aortic aneurysm). An error metric was defined to quantify the quality of the VDM results. Furthermore, a k-means clustering technique was used to categorize the subjects into three clusters based on ascending aortic motion. Optimal penalty parameters were identified for each of the three clusters. The results indicated that patient clusters with smaller aortic root motion required larger rigidity penalty values. The calibrated parameters successively reduced errors in 3D displacement and multi-axial stretch compared to un-optimized VDM predictions, enhancing the accuracy of capturing aortic deformation from dynamic images. Among the different aortic regions, the ascending thoracic aorta exhibits the largest error reduction.

The Poisson-Nernst-Planck (PNP) system is an accurate model of electrodiffusion of ionic species. It is commonly used in situations where nanoscale resolution is required, for instance close to ion channels in the membranes of biological cells. The inherent stiffness of the equations has made them challenging to solve and has limited the applicability of the system. In particular, the time step required for stable solutions has typically needed to be very short (nanoseconds), which makes simulations on the time scale of an action potential (milliseconds) difficult. Recently, it has been observed that avoiding operator splitting and instead solving the concentration equations and the electrostatic equation in a coupled manner relaxes the time-step limitation considerably. However, no theoretical explanation of this observation has been provided. Here, we aim to explain why the coupled scheme allows much larger time steps. We illustrate the mechanism by considering special cases that define necessary, but not sufficient, conditions for stability. We also show that these conditions remain relevant for the fully coupled PNP model in 3D.

Uncertainty quantification (UQ) in computational heart models is essential for reliable cardiac digital twins (DTs) in personalized medicine, yet remains challenging. Traditional Monte Carlo and stochastic Galerkin methods often become impractical in the high-dimensional, nonlinear state variable and parameter spaces of cardiac electrophysiology and mechanics. This article introduces a framework for learning a joint probability density over cardiac observables and model parameters, enabling the characterization of statistical dependencies across a large number of variables in patient-specific cardiac DTs. By sampling from this density and conditioning on available data, useful predictive distributions can be constructed, allowing uncertainty to be propagated through the model and quantified in terms of variability. Conditional regression can then be performed directly on this learned density, enabling systematic exploration of interdependencies among observables for both predictive inference and model design. The statistical methodology adopts a geometry-aware generative learning framework, recently introduced by the authors, that decouples the learning of data geometry from sampling. First it identifies a low-dimensional latent representation that captures the intrinsic structure of the data and its multiscale geometric features. A stochastic differential equation is then formulated directly in the low-dimensional latent space to generate samples efficiently; these are subsequently mapped back to the high-dimensional space of cardiac states and parameters through a smooth lifting operator. We demonstrate the approach on a ventricular arrhythmia prediction benchmark, where the learned joint probability density enables the construction of predictive distributions over key parameters (e.g., conductivities, fibrosis patterns) through sampling and conditioning. This enables uncertainty to be propagated and quantified through sampling and conditioning on the learned joint density, with substantially fewer model evaluations than conventional UQ methods.

A major challenge in neuroscience is to predict how currents in nanodomains affect voltage and ionic concentrations. Cable and Rall theory provide analytic current-voltage relations by neglecting concentration gradients, and the impact of concentration gradients is usually studied numerically with the Poisson-Nernst-Planck (PNP) model. A precise quantitative understanding of the combined dynamics remains limited because analytic current-voltage-concentration relations are missing. In this work we derive such relations using a novel approach based on cross-diffusion equations. For narrow cylindrical domains, we derive time-dependent and steady-state expressions that explicitly show how currents affect voltage and ionic concentrations. We find that the influx of only one ion can significantly change the concentrations of all the other ions even if no channels for these ions are present. After a current injection we compute a biphasic voltage transient where the small-time asymptotic corresponds to the steady-state solution of the cable equation. We show that the accuracy of cable theory prediction for the voltage depends on how the current is distributed among the various ions. Finally, we develop an iterative method to accurately compute steady-state profiles for voltage and concentrations using first-order results by subdividing a cylinder into small segments.

Hemodynamic forces play a key role in early cardiac morphogenesis, yet many computational studies assume Newtonian blood behavior. Here, we evaluate the impact of nonNewtonian shearthinning rheology on flow patterns, pressure distributions, and wall shear stress (WSS) during cardiac looping using idealized threedimensional models of the embryonic heart tube. Five geometries representing progressive looping stages, from a linear tube to an Sshaped configuration with ventricular ballooning, were analyzed under pulsatile flow using both Newtonian and powerlaw viscosity models. Across all stages, Reynolds numbers (Re {approx} 1-7) and Womersley numbers (Wo {approx} 0.3) indicated laminar, quasisteady flow consistent with embryonic conditions. Incorporating shearthinning rheology produced substantial deviations from Newtonian predictions, with peak systolic WSS differing by up to [~]40% and pressure drops by up to [~]20%. These effects were most pronounced in regions of increased curvature and geometric complexity. These findings demonstrate that nonNewtonian rheology significantly influences predicted hemodynamic environments during cardiac looping and should be incorporated into computational models aimed at understanding mechanobiological regulation of early heart development.

The aorta, normally resilient to hemodynamic stresses, becomes vulnerable to structural failure due to diverse conditions that weaken the wall. We injected fluid into excised specimens of human ascending aorta with pressure monitoring to quantify the impact of clinical and histological factors on mural damage. Two modes of medial injury emerged with distinct pressure tracings. Extravasation was characterized by diffuse infiltration of fluid with widespread damage of smooth muscle cells and collagen fibers but limited separation of elastic lamellae. By contrast, delamination was characterized by marked separation of elastic lamellae along a single plane with damage to cells and fibrillar matrix restricted to adjacent laminae. Aging, aortic dilatation, and family history associated with lower pressures causing delamination, whereas a diagnosis of hypertension associated with higher pressures suggesting resilience to dissection. Collagen fraction adjacent to delamination correlated with higher pressures as did decreased smooth muscle cell density and increased glycosaminoglycan fraction, although several clinical and histological variables were interrelated. Protein cross-linking strengthened and enzymatic digestion of collagen weakened the aortic wall, while acute cell lysis with detergent had no effect. We conclude that increased functional medial collagen has an adaptive protective role in aortic remodeling rather than signifying medial degeneration.

Arterial thrombosis is initiated when mechanical forces in flowing blood exceed the activation thresholds of platelets and von Willebrand factor (vWF). Despite extensive experimental characterization of shear-induced platelet aggregation, a unified theoretical framework that maps hemodynamic forcing onto clot nucleation is lacking. Here we present Force-Gated Thrombosis (FGT), a non-equilibrium mechanical theory that treats thrombus formation as a continuous phase transition driven by an effective mechanical forcing {Sigma} ={sigma} + |{nabla}{sigma}| + {beta}{varepsilon}, which combines local wall shear stress{sigma} , shear gradient |{nabla}{sigma}|, and extensional strain rate{varepsilon} . We introduce a dimensionless Thrombosis Number {Theta} = ({Sigma}/{Sigma}c)(P/P0)m(C/C0)n, which incorporates platelet concentration P and coagulation factor concentration C, and governs the transition between stable flow ({Theta} < 1) and active clot growth ({Theta} > 1). The thrombus density is represented by a scalar order parameter{varphi} whose dynamics follow a Ginzburg- Landau free energy functional. For a simplified stenosed artery we derive an analytic closed-form thrombosis onset criterion and a critical flow rate [Formula], where{delta} is stenosis severity. Linear stability analysis shows that perturbations grow at rate{omega} (k) = {Lambda}({Theta}) - D{varphi}k2, becoming unstable when {Theta} > 1. Near threshold the clot volume fraction scales as{varphi} [~] ({Theta} - 1)1/2, a mean-field critical exponent consistent with Ginzburg- Landau theory. Systematic comparison with fifteen published experimental and computational datasets spanning shear rates from 100 to 15,000 s-1 confirms that FGT correctly predicts the existence, location, and approximate severity of pathological thrombus formation across diverse vascular geometries. The theory provides a quantitative bridge between single-molecule mechanobiology and macroscale clinical thrombosis, and yields experimentally testable predictions distinguishing FGT from purely biochemical models.

Vertex models are widely used within the field of developmental biology to study tissue morphogenesis. These models are well-suited for modeling deformation at the cellular level where movement is driven by local forces. However, understanding how these microscopic movements coordinate to yield macroscopic phenomena such as the shapes of entire tissues remains a challenge. Here we study a top-down approach using differentiable programming on a simplified vertex model of a laminar tissue, and investigate whether the attributes of individual cells can be tuned to make the mesh as a whole acquire a predefined shape. We let the mesh evolve according to simple rules defined by the input to each polygon, and evaluate the resulting shape against a target boundary. Additionally, we show how the high degeneracy of the output can be reduced by constraining the polygon distributions: first, by adding simple penalties on tissue-wide attributes; and second, by dividing the tissue into regions, within which we bias the attributes toward characteristic values. Our study shows how a simple vertex model can be combined with differentiable programming to model developing tissues, and provides insight into the way individual cells must coordinate to yield macroscopic phenomena such as pre-programmed shapes.

Quantifying the diagnostic dispersion of inferred parameter distributions is a challenge in uncertainty-aware modeling. Scalar summaries such as credible interval width are topology-blind; fundamentally different posterior morphologies can yield identical scores, obscuring whether a parameter is precisely estimated or constrained to a range. We propose a Composite Certainty Framework that addresses this metric degeneracy by aggregating five complementary uncertainty metrics including interquartile range, standard deviation, full width at half maximum, Shannon entropy, and mass width. These metrics are aggregated through non-parametric Borda rank voting into a single, unitless consensus certainty score. Applied to a simulation-based inference pipeline for a finite-element model of the human middle ear tuned to cadaveric acoustic measurements, the framework reveals parameter-specific identifiability profiles invisible to any individual metric. It produces two actionable clinical thresholds: (1) the maximum tolerable measurement noise for reliable parameter recovery, and (2) the minimum simulation budget for posterior convergence. We demonstrated that no single metric captures all aspects of posterior dispersion, as spread-based metrics and entropy diverge systematically for clinically critical parameters, whereas their aggregation produces a consensus reflecting genuine diagnostic certainty. The framework is generalizable to any model-based diagnostic pipeline where posterior distribution not merely its coverage, but determines clinical certainty.

Background and ObjectiveThe foot-ankle complex is a highly articulated and mechanically constrained system, often simplified as a chain of few rigid segments, neglecting many bone-to-bone motions and raising questions about the accurate representation of interaction with ground. This study proposes a new reduced-order multibody formulation that captures intrinsic kinematic constraints of the foot through motion synergies. MethodsBones kinematic coupling, or motion synergies, were experimentally derived from weight-bearing CT scans using principal component analysis. These couplings were embedded in a synergy-based multibody kinematic optimization framework describing the foot-ankle with five degrees of freedom: ankle flexion; foot adduction, pronation, and arching; and toe flexion. Model accuracy was evaluated against bone-level experimental kinematics. The model was applied to gait data from healthy, flat, and diabetic feet and compared with a standard multi-segment foot model, assessing robustness by progressively reducing the number of skin markers. ResultsAverage errors were about 1{degrees} and 0.5 mm when using subject-specific synergies and below 7{degrees} and 4 mm when scaling the generic model, matching or exceeding the accuracy of existing models. Reliable reconstruction was obtained using only four foot markers. In clinical gait analysis, the model showed superior discrimination between populations and enabled assessment of transverse arch deformation, not accessible with conventional models. ConclusionThe proposed synergy-based model provides an accurate, low-complexity framework for reconstructing bone-level foot and ankle kinematics, substantially simplifying gait analysis while improving biomechanical interpretability. This framework supports future integration with dynamic models aimed at studying load transmission in the foot.

Curvature evolution on a deforming surface is governed by the full change in the surface metric, but on biological surfaces captured by serial three-dimensional imaging, only the local area change is observable. The loss of the shear component leaves prediction of curvature evolution underdetermined from imaging alone. On the thoracic aorta, where curvature change marks disease progression, we derive a closed-form equation that predicts the change in integrated Gaussian curvature from the area dilation and initial geometry. The equation combines a conformal term in the area dilation with a leading anisotropy correction from the initial geometry. These two analytic levels, augmented by multi-scale spatial features at neighboring regions and a graph neural network trained on residuals, form a four-level nested predictor. On a synthetic aortic geometry under prescribed isotropic expansion, the equation recovers the analytic coefficient exactly. Across a continuum from pure expansion to pure shear, it holds R2 [&ge;] 0.71. On 236 paired thoracic aortic surfaces spanning dissection, aneurysm, traumatic injury, and non-pathologic controls, the equation recovers within-surface curvature change patterns with per-patient median Pearson [Formula] and pooled R2 = +0.238 [+0.225, +0.250], matching the graph neural network on the same inputs. The residual is a direct measurement of how far the observed growth field departs from conformality. HighlightsO_LIClosed-form equation predicts aortic curvature change from paired computed tomography scans. C_LIO_LIRecovers analytic predictions exactly on synthetic aortic geometries. C_LIO_LIAnisotropy proxy holds R2 [&ge;] 0.71 from pure expansion to pure shear. C_LIO_LICoefficients tie to geometric mechanisms ensuring interpretability. C_LIO_LIAnisotropy term, computable from one CT, is twice as large on diseased aortas. C_LI

Biomechanical simulations of complex hand motions remain scarce, due to challenges that span computation and data acquisition. Using a computer vision-based motion capture approach, a 23-degree of freedom musculoskeletal model, and direct collocation optimization, we performed muscle-driven simulations to track hand kinematics from 7 participants performing American Sign Language gestures. While proximal joints were tracked accurately, interphalangeal joint tracking was significantly worse, with a consistent flexion bias. Modifications to finger extensor muscle paths that incorporated the dual-inserting nature of the extensors improved accuracy, suggesting better representation of extensor force distribution across distal joints may be necessary for accurate hand simulations.

The human myometrium undergoes a dramatic transformation during pregnancy, shifting from quiescence to highly synchronised contractility. Understanding this transition is crucial for addressing pathologies such as preterm labour and dystocia (ineffective labour). We present a multi-scale Functional Tissue Unit (FTU) model allowing us to investigate how tissue-level excitability emerges from single-cell electrophysiology. We propose a heterogeneity-driven selection mechanism, wherein a sub-population of cells with high intrinsic excitability dynamically emerges as pace-makers. This active process complements passive depolarisation by interstitial cells, allowing spontaneous excitation to arise without a fixed anatomical pacemaker. Stochastic simulations produced an average burst frequency of 0.047 Hz ({approx}2.8 bursts per minute), closely consistent with clinical measurements of 2-3 contractions per minute during active labour, and demonstrated that this function is robust to spatial topological changes. Furthermore, implementation of inflammation-induced remodelling simulations successfully linked molecular-level changes to a preterm labour phenotype. This model provides a platform for investigating uterine contractility and serves as a component for future whole-organ Physiome models.

Mechanical properties of epithelial tissues play essential roles in morphogenesis and physiological function. In this study, we analytically derived the in-plane bulk modulus, shear modulus, and Poissons ratio of a three-dimensional cell vertex model of epithelial monolayers. We showed that the model can robustly reproduce a near-zero in-plane Poissons ratio, a mechanical feature reported in cultured epithelial tissues. Numerical simulations further confirmed that the theoretically predicted Poissons ratio accurately describes the response of the model under finite, biologically relevant strains. In addition, the model exhibits not only morphological bistability between squamous-like and columnar-like states, but also mechanical bistability characterized by distinct elastic responses. Together, these results provide a minimal three-dimensional framework that links cell-scale mechanical interactions and epithelial morphology to tissue-scale elastic properties.

Patellar osteochondral allograft (OCA) transplantation is widely used to treat large full-thickness cartilage defects, yet long-term failure and reoperation rates remain high. Although surface congruity and osseous integration are emphasized clinically, cartilage thickness and mechanical compatibility between donor and recipient are not considered. Our previous work suggests that cartilage thickness mismatch can amplify local deformation at the graft boundary, potentially compromising graft longevity. This study investigates how combined mismatches in cartilage thickness and mechanical properties influence the local strain environment at the patellar OCA interface. Simplified two-dimensional axisymmetric finite element models of patellar OCA repair were developed in ABAQUS. Donor-to-recipient cartilage thickness ratios ranging from 0.33 to 3.25 were evaluated together with donor-recipient Youngs modulus mismatches (2.5-7.0 MPa). Cartilage was modeled using homogeneous linear elastic and functionally graded material formulations to account for depth-dependent stiffness. A compressive pressure of 1.0 MPa was applied to represent patellofemoral joint loading, and peak compressive and shear strains were quantified at the graft boundary. Cartilage thickness mismatch produced localized high-strain regions (HSR) of compressive and shear strain at the donor-recipient interface that were absent in thickness-matched constructs. Strain amplification increased with both thickness and mechanical property mismatch. Compressive strain exhibited directional asymmetry, with donor-side-thicker configurations producing greater amplification than recipient-side-thicker configurations. Incorporating depth-dependent cartilage stiffness reduced peak strain magnitudes but did not eliminate mismatch-driven strain amplification. These findings demonstrate that cartilage thickness and mechanical disparity can create HSR at the patellar OCA graft boundary that may predispose grafts to impaired integration and long-term failure.

Understanding how biological shape and movement interact with surrounding fluids represents a fundamental challenge at the intersection of biology, physics, and engineering. Fish locomotion exemplifies this challenge: body morphology and swimming kinematics together determine the hydrodynamic forces and flow structures that enable efficient propulsion and maneuverability. Whereas biologists have long sought to connect morphological variation to swimming performance, traditional morphometric approaches provide limited insight into the fluid mechanical consequences of shape differences. Similarly, although computational fluid dynamics can reveal detailed flow physics, simulating hydrodynamics across diverse and dynamic morphologies remains prohibitively expensive for systematic investigation. To bridge this gap, we introduce a data-driven framework that connects fish body shape dynamics to hydro-dynamic performance through compact morphospace parameterization and reduced-order modeling. Using CFD simulations of 15 fish species from the Digital Life Project database (www.digitallife3d.org/3d-model), we generate hydrodynamic datasets capturing the shape-flow relationship. Principal Component Analysis (PCA) extracts four dominant shape parameters from dorsal body profiles, which are then integrated into an Inverse-Design with Dynamic Mode Decomposition (ID-DMD) framework to model the resulting fluid dynamics. The resulting modal analysis suggests that locomotion strategies emerge from specific shape-flow interactions. We further demonstrate the frameworks utility through single- and multi-objective shape optimization, showing how it enables efficient exploration of the morphology-hydrodynamics relationship. This approach offers a novel analysis and design tool for understanding how biological form and motion interact with fluid mechanics, with applications ranging from bio-inspired vehicle development to evolutionary biomechanics.

Idiopathic pulmonary fibrosis (IPF) is a progressive and ultimately fatal disease of aging, driven by dysregulated fibroblast activation and accompanied by collagen accumulation in the lung interstitium, resulting in tissue stiffening. While the accumulation of senescent cells has been increasingly implicated in IPF pathogenesis, understanding the reciprocal dynamics of senescent fibroblast levels and evolving tissue mechanics is difficult to achieve with experimental approaches alone. To address this limitation, we developed an agent-based model (ABM) of fibroblast activation in the lung that couples cell behavior to the dynamic mechanical changes accompanying fibrosis. This model was parameterized entirely from experimental data in young mice to enable robust validation and then adapted to fit aged mouse biology for additional validation. Both young and aged models accurately reflected changes in collagen accumulation and stiffness burden of experimental systems. We then incorporated senescent cell behavior into the aged model to investigate how senescent cell burden influences fibrosis progression and how cell-cell interactions drive senescent cell accumulation. These simulations identified a unique role for juxtacrine-mediated contact between non-senescent and senescent fibroblasts in expanding the total senescent cell burden. Our ABM also revealed that the timing of immune-mediated senescent cell clearance critically regulates fibrotic outcomes. Together, this ABM provides useful insights into how the interrelated dynamics of tissue mechanics and senescent fibroblasts drive fibrosis progression.

Take-off is a fast and energy-efficient strategy for bipedal animals, such as birds, to achieve rapid movement; however, how muscle physiology scales to govern this universal behavior remains unresolved. Research in other species physiologies is not readily applicable. As a result, important questions, whether theropod dinosaurs such as Tyrannosaurus rex were capable of jumping, remain unanswered. In this article, we coupled Lagrangian dynamics with Hills muscle equations and developed new experimental methods to quantify joint rotational stiffness and damping, thereby enabling a systematic description of lower-limb mechanics. The approach establishes a novel kinetic framework that links muscle contractile properties to lower-limb performance without invoking control optimization. Animal observations and tabletop mechanisms validate the framework. The mechanics model reveals that the take-off time of about 0.1 s across body masses of 0.003 to 90 kg is achievable, as heavier birds generate proportionally higher reaction forces. Additionally, Tyrannosaurus rex should be capable of jumping, based on the available physiology data. Beyond evolutionary insights, our framework provides a new methodology for analyzing the mechanical properties of biological joints and informing the design of scalable bio-inspired robots.